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ABSTRACT 


A fast multigrid solver for the steady, incompressible Navier-Stokes equations is presented. The 
multigrid solver is based upon a factorizable discrete scheme for the velocity-pressure form of the 
Navier-Stokes equations. This scheme correctly distinguishes between the advection-diffusion and 
elliptic parts of the operator, allowing efficient smoothers to be constructed. To evaluate the 
multigrid algorithm, solutions are computed for flow over a flat plate, parabola, and a Karnran- 
Trefftz airfoil. Both nonlifting and lifting airfoil flows are considered, with a Reynolds number 
range of 200 to 800. Convergence and accuracy of the algorithm are discussed. Using Gauss-Seidel 
line relaxation in alternating directions, multigrid convergence behavior approaching that of 0(N ) 
methods is achieved. The computational efficiency of the numerical scheme is compared with that 
of Runge-Kutta and implicit upwind based multigrid methods. 
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1 INTRODUCTION 


One of the critical needs in computational fluid dynamics is faster flow solvers. Multigrid is a 
well known method of convergence acceleration that is widely used in Euler and Reynolds-averaged 
Navier-Stokes codes. These applications of multigrid generally are based on the unsteady equations 
using some temporal integrator as the smoother, combined with a full-approximation scheme (FAS) 
multigrid iteration. A common approach is one originally proposed by Jameson [1]. Starting with 
the unsteady equations, a finite-volume spatial discretization with explicit artificial viscosity is 
combined with a Runge-Kutta (R-K) time integration as a smoother. An alternative approach 
[2, 3, 4] is to use upwind-differencing and implicit time integration as the smoother. However, these 
approaches have resulted in poor multigrid efficiency. When applied to high Reynolds number flows 
over complex geometries, convergence rates are often worse than 0.99. There is clearly a need to 
develop substantially more efficient multigrid solvers. 

According to Brandt [7], one of the major obstacles to achieving better multigrid performance for 
advection dominated flows is that the coarse grid provides only a fraction of the needed correction for 
some smooth error components. This obstacle can be removed by designing a solver that effectively 
distinguishes between the advection and elliptic (hyperbolic if Mach number exceeds one) factors 
of the system and treats each one appropriately. For instance, advection can be treated by space 
marching, while elliptic factors can be treated by multigrid. The efficiency of such an algorithm will 
be essentially identical to that of the solver for the elliptic factor only, and thereby attain so-called 
“textbook” multigrid efficiency. Such efficiency can also be viewed as the efficiency that can be 
attained when solving a pure elliptic problem (e.g., Poisson equation). 

One essential element for achieving textbook multigrid efficiency is an effective multigrid scheme. 
The performance of a multigrid scheme strongly depends upon two components: 1) relaxation 
(smoother), 2) coarse-grid operator. A relaxation process is devised such that it effectively smooths, 
as determined by its smoothing factor, the high-frequency components of the solution error, and a 
coarse-grid operator is designed so that it efficiently reduces the low-frequency error components. 
Generally, a smoothing factor of 0.5 is considered good. Such a scheme is said to have a good 
measure of h-ellipticity. Often, the amenability of a discrete operator to a multigrid algorithm is 
characterized by the h-ellipticity of the operator. For the purpose of subsequent discussion the 
definiton of h-ellipticity is given here. The measure of h-ellipticity is defined in the following way 
(see discussion in Trottenberg, Oosterlee, and Schuller [5]). Let L h represent the discrete operator 
matrix for the system of equations on a discrete domain with uniform spacing h. Also, let E \ be 
the measure of h-ellipticity. Then, 

min ( | det t h (9 x ,9 y )\ : f < \0 X \, \9 y \ < 7r) 

E h { Lh) = V f , 

max 1 1 det L h (9 x ,9 y )\ : 0 < \0 X \, \9 y \ < 7rj 

where L/, is the Fourier symbol of the operator L h, and 9 X , 9 y are proportional to the wave numbers 
associated with the x and y coordinate directions, respectively. According to this definition, good h- 
ellipticity means that Eh must be bounded away from zero by a suitable constant. As a reference, a 
standard five-point discretization of the Laplacian has an Eh of 0.25. To achieve textbook multigrid 
efficiency, the discrete operators of the multigrid algorithm must have good h-ellipticity. 

Another essential element for attaining textbook multigrid efficiency is factorizability. A scheme 
is defined to be factorizable if the determinant of the discrete operator matrix can be represented 
as the product of separate scalar factors (e.g., advection-diffusion, elliptic). With a factorizable 
scheme the solution to the system is reduced to the solution of the separate scalar factors. The key 
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aspect of the factorization is that only those terms that contribute strongly to the high-frequency 
components of the error must be considered when constructing the relaxation method. These terms 
are called the “principal part” of the operator. The approach advocated by Brandt is to introduce 
a transformation of the dependent variables to a set of auxiliary variables, which he calls “ghost” 
variables, such that the principal part of the operator is triangular. The equation residuals produce 
increments in the ghost variables, which are distributed to the primary variables by means of the 
transformation (distribution) matrix. This is called “distributive relaxation” (for further discussion 
see Brandt [7], Hackbush [6], Trottenberg et al. [5]). An alternative to distributive relaxation is to 
find a projection matrix for the original system such that the principal part of the resulting discrete 
system is in triangular form, but in terms of the orignal set of dependent variables. This is the 
approach used here. 

Relatively little research has been done in the area of multigrid algorithms for the Navier- 
Stokes equations based upon factorizable discrete schemes. Using a scheme based on distributive 
relaxation, Brandt and Yavneh [8] have demonstrated textbook multigrid efficiency for the incom- 
pressible Navier-Stokes equations for entering flow. Their results are for a simple geometry and a 
Cartesian grid, using a staggered-grid discretization (i.e., velocity components and pressure stored 
at different locations) of the equations. 

Recently, Thomas, Diskin, and Brandt [9] achieved textbook multigrid efficiency for high Reynolds 
number incompressible wake and boundary layer flows associated with a flat plate. Their scheme 
uses a staggered grid approach with distributed relaxation and defect correction. With the dis- 
tributive relaxation the system of equations was decomposed (i.e., factorized) everywhere, except 
near boundaries where they remained coupled. In all calculations Cartesian grids were used. 

Sidilkover and Asher[10] introduced a fast multigrid solver for the incompressible Navier-Stokes 
equations that does not require a staggered grid arrangement of flow variables. With a pressure 
Poisson formulation as a foundation, a factorizable discrete scheme is derived. In this work also 
only simple model problems were solved using uniform Cartesian grids. 

In this paper a generalization of the approach of Sidilkover and Ascher[10], which is a continua- 
tion of the work of Roberts and Swanson [11] and Swanson [12], is used. A conventional vertex-based 
finite- volume or finite-difference discretization of the primitive variables is used, avoiding the need 
for staggered grids. This simplifies the restriction and prolongation operations, because the same 
operator can be used for all variables. A projection operator is applied to the incompressible 
Navier-Stokes equations in velocity-divergence form to obtain a factorizable discrete scheme with 
decoupling for the velocity-pressure form of the Navier-Stokes equations, except near the bound- 
aries of the domain being considered. For the pressure equation a numerical a boundary condition 
based on the momentum equations is introduced to enforce continuity at the inflow and solid surface 
boundaries. Gauss-Seidel relaxation is used as a smoother of the Poisson equation for the pressure. 
The momentum equations are relaxed by space marching. Since the elliptic and advection-diffusion 
parts of the system are decoupled over most of the flow field, ideal multigrid efficiency is possible. 
With the present scheme line relaxation is used to treat the anisotropy of mesh aspect ratio that 
generally occurs in resolving viscous flows. 

The primary purpose of this paper is to evaluate a modified version of the scheme introduced by 
Swanson [12], herein called INCOMP scheme. In Section 2 the governing equations are presented. 
An iterative operator matrix is determined by applying a projection matrix. Discretization of 
the flow equations and the associated boundary conditions are discussed in Sections 3 and 4. 
Section 5 describes the multigrid algorithm for solving the equations. Then, in Section 6 there is 
discussion about computation of skin-friction drag. Surface skin friction is as an important quantity 
in assessing a viscous flow solver. Next, computed results for three types of flows are shown and 
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discussed in Section 7. Flow solutions were calculated for the following geometries: 1) flat plate, 
2) parabola, 3) airfoil. Both accuracy and convergence behavior are considered in evaluating the 
present scheme. The computational work required by the INCOMP scheme is compared to that 
needed by R-K and implicit upwind multigrid solvers. Concluding remarks are made in Section 8. 

2 MATHEMATICAL FORMULATION 

The boundary-value problem for the two-dimensional incompressible Navier-Stokes equations can 
be written as 


(v • v) v + Vp - i/Av = f , x g n 

V-V = 0, xeH (1) 

B(v,p) = o, xedn 

where V = V(x, y) = ( u , v) T is the fluid velocity, p = p(x, y) is the pressure, and A is the Laplacian 
operator. The density is taken to be one, and the coefficient v is the reciprocal of the Reynolds 
number, v = 1/Re. The domain II is in 3ft 2 , and the boundary of the domain is dfl. There are two 
boundary conditions, and these are represented by £>(V, p). In Section 4 the boundary conditions 
will be discussed in detail. Throughout this paper the forcing function f = f (x,y) = 0. The system 
of Eq. (1) is sometimes called the velocity-divergence form of the equations. 

The velocity-divergence formulation of Eq. (1) can be used with a staggered-grid discretiza- 
tion. However, such a formulation is not appropriate for a collocated-grid discretization when the 
continuity equation V • V and pressure gradient Vp are approximated with central differencing. 
Under these circumstances spurious modes in the pressure can be produced (i.e. , there is a loss of 
h-ellipticity) . 

There are alternative forms for the system of Eq. (1) that are amenable to collocated-grid 
discretization and to a rapidly convergent numerical soluion procedure, such as one based on a 
muligrid scheme. One of these forms can be derived by taking the divergence of the momentum 
equation and applying V- V = 0. The resulting system is is called the velocity-pressure (or pressure 
Poisson) form of the equations, and the associated boundary- value problem is defined by 


(V • V) V + Vp - AV = f , x g n 

A p + Vu • \ x + Vu • \ y = V • f , x € fl 

B(V,p) = 0, xGffl 

V • V = 0, xGffi, 


(2) 


Now, the pressure equation plays the role of the continuity equation in Q. Since this equation is 
second order, the system requires an additional boundary condition. By choosing the additional 
boundary condition to be V-V = 0, one can prove that the velocity-pressure and velocity-divergence 
formulations are equivalent (see Trottenberg et al. [5]). 

The velocity- pressure form of the equations given in Eq. (2) can be derived other ways that 
can clearly suggest how to construct a fast solution algorithm. One way is to start by defining an 
advection-diffusion operator 

Q V = Q- uA, (3) 

where Q = ud x + vd y , and d x , d y are the partial differentiation operators. Equation (1) can then 
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be written as 


(Qv 0 d x \ (u\ 

Lq = 0 Q u d y I u 1 = 0. (4) 

\8 X Oy 0 / \pj 

Next define the operator Q u according to 

QAf ) = Q*(f) + a (uf) = —d x {uf) - d y ( v f ) + A(i //), (5) 

where Q* is the adjoint of Q, and construct a projection operator P as 


(l 0 0\ 

P = o 1 0 . ( 6 ) 

\d x dy Q u ) 

After applying the projection operator to the velocity-divergence form of Eq. (1), which is equivalent 
to taking the divergence of the momentum equation and applying V • V = 0, one obtains 


(Qu 0 d x \ (u\ 

PLq = 0 Q u d y 1 u I - b = 0, 

\0 0 Aj \pj 

or 

Lq = b 

where b = (0 0 63 ) T and 

&3 = d x ( v Uy - UVy) - dy(uv x - VU X ) 


( 7 ) 

( 8 ) 
(9) 


The 3x3 upper triangular operator matrix of Eq. (7) corresponds to the operator obtained by 
linearizing the velocity-pressure form of the system and retaining only principal terms. The princi- 
pal terms are those that dominate with respect to high-frequency error components. The column 
vector b is comprised only of subprincipal terms over most of the domain, and these terms are not 
important for the purpose of constructing a relaxation scheme. 


3 DISCRETIZATION 


The first step in approximating L is to discretize the Navier-Stokes equations given in Eq. (4). 
Consider a typical grid vertex (i,j) as shown in Fig. 1. The equations are discretized with a 
combination of finite-difference and finite- volume techniques. For the momentum equations second- 
order accurate upwind differencing is used for the advection part of the operator Q u in Eq. (3). 
Central differencing is used for the physical diffusion contribution to Q u and the pressure gradient 
term. Discretization of the pressure equation in Eq. (7) is also based on central differencing. 

Consider the discretization of the advection term Qu in the x'-rnornentum equation. On a 
Cartesian grid the second-order upwind-difference discretization of the contributions to Qu is 


ud x u\ij — ^j^s u \uij\(3 4 E x u + E x U )ui : j, 
vd y u\ij = I'D ,j \ (3 v + Ey l )uij, 


(10) 
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Figure 1: Primary and dual cells on quadrilateral grid. 


where the shift operator E is defined by 


UiJ — 'U'i— 


i—s u m,j 5 


E~ 


/ U j i,j — 'U'i.i— 


i,J—s v rm 


(ii) 


with s u = sign(uij ) and s v = sign(vi } j). Here, the superscript h denotes the discrete approximation 
to the corresponding differential operator. For arbitrary curvilinear coordinates (£, rj) the advection 
operator Q is defined by 

Qu = (ud x + vd y )u = (Ud^ + Vd v )u ( 12 ) 

where U = £ x u+(,yV and V = rj x u+r] y v are contravariant velocity components. The transformation 
from Cartesian coordinates (x,y) to generalized coordinates (£, 77 ), which correspond to the (i,j) 
indices of the grid, is given by 


(f 7h ) = J " 1 ( Vr > , (13) 

\iy VyJ V X V x £ J 

and the transformation Jacobian J = x^y v — x v y^. The terms JJd^u and Vd v u are differenced 
according to Eq. (10) on a uniform computational space. 

To approximate the physical viscous terms in the momentum equations we first compute the 
gradients of the velocity components at the edges of the dual grid cell, delineated by dashed lines in 
Fig. 1. The evaluation points are shown as solid squares. For example, the gradients of the velocity 
component u at the point (i — 1 / 2 , j) of the dual grid cell are 


u x — T UrjVx: y T U y 1)y, 


(14) 


The derivatives u % and u v are approximated at the face center by 


u (\i-l/2,j ~ u h3 u i-l ,ji 

1 1 , , 

— /j_ (^*,1+1 "F 1,1+1 1 V'i—lJ—l) j 


(15) 


with similar expressions at the points ( i + 1/2, j), ( i,j + 1/2), and ( i,j — 1/2). The gradients of 
pressure (p) are calculated the same way. The grid metric terms £ x , £ y , rj x , and r] y are also evaluated 
on the dual grid face centers. 
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With the velocity gradients known at the faces of the dual grid cell, one can easily approximate 
the Laplacian in the viscous terms Au and Au. Applying Green’s theorem to the dual grid cell, 


An = — — * (■ u x dy — Uydx ), Av = — — ® (v x dy — v y dx), 

dAi.j dAi.-i 


(16) 


where Ay is the area of the dual cell at the vertex (i. j). 

The projection operator P is applied to the discrete equations to obtain the residual for the 
pressure Poisson equation of Eq. (7). In deriving Eq. (7) we assumed that the functions u(x,y) 
and v(x,y), as well as their spatial derivatives, are continuous in 17 and on 917. Thus, the order 
of differentiation does not matter, and the viscous terms are eliminated in the pressure equation. 
Here, we assume that the discrete operators exhibit this property. Thus we are assuming that the 
differential equalities 

dy dy d x = d x dy dy, d x d x dy = dy d x d x . (17) 

have discrete analogs. Letting ( r p )ij be the pressure equation residual at vertex (i,j), the applica- 
tion of P can be written in integral form, 


(' r p)i,j ~ 


dA, 


Q h u + d^p 


— u 


d h 


u 


dS 


dy - ( 


Q h v + d h vP 


— V 



dx. (18) 


The discretization of the pressure equation is done by first discretizing the boxed terms of Eq. (18) on 
the edges of the dual grid cell. The integral of Eq. (18) is then evaluated to get the pressure equation 
residual at the vertex ( i , j), taking the boxed terms to be constant over each face of the dual grid cell. 
For a uniform Cartesian grid, the principal part of the resulting discretization is a conventional five- 
point approximation to the Laplacian operating on the pressure. At inflow and surface boundaries 
additional terms appear in the residual due to the implementation of the divergence-free boundary 
condition. These terms will be defined in Section 4 on boundary conditions. 


4 BOUNDARY CONDITIONS 

For a given physical problem let the domain of interest and its boundary be denoted by 17 and <917, 
respectively. Let the boundary be partitioned according to 

d£l = Ti n + T ou t + T sur f, (19) 

where E; n is the inflow part of the boundary, r o . U f is the outflow part, and T sur f is the solid surface 
part. For each T the velocity-pressure formulation of the incompressible N-S equations requires 
three boundary conditions. The boundary conditions for the velocity- divergence formulation are 
augmented by an additional boundary condition associated with the pressure equation. 

There has been considerable discussion in the literature (e.g., [14] - [18]) regarding the proper 
boundary condition for the pressure equation. According to Gresho and Sani [14] the normal com- 
ponent of the momentum equation is the appropriate boundary condition. However, as indicated 
by Henshaw [15], this particular boundary condition allows a nonzero divergence on the boundary, 
and thus, allows a nonzero divergence in the domain 17. Furthermore, Strikwerda [17] points out 
that this boundary condition provides no new information, resulting in an underdetermined system. 
The viewpoint adopted here is that V • V = 0 is the proper boundary condition. With V • V = 0 on 
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the boundary dfl and sufficiently smooth solutions, the velocity-divergence and velocity-pressure 
formulations are equivalent. This leads to a well-posed problem. The momentum equation provides 
a suitable vehicle for implementing this boundary condition in the discrete problem. 

The boundary conditions along the boundary d£l are specified according to 


u(x, y) = U in (x , y ), v(x, y) = v in (x, y) 

on 

r. 

L in 

(20) 

V • t = Vt(x, y), p(x,y) =Pout(x,y) 

on 

r out 

(21) 

u(x,y) = 0, v(x,y) = 0 

on 

r surf 

(22) 

< 

< 

II 

0 

on 

dfl 

(23) 


where t is a unit tangent vector and the subscript t means tangential component. 

The normal component of the momentum equation can be used to implement the divergence- 
free boundary condition. Taking the inner product of the unit normal to the boundary and the 
vector momentum equation, this equation is 


dp 

dn 


EE 'W'kVjdjVk V EE n kdjdjVk, 

j = 1 k = 1 j = 1 k = 1 


(24) 


where V\ = u, V 2 = v, and the Cartesian components for the inward pointing unit normal to a 
boundary are n\ = n x and 112 = n y . In the discrete problem the normal pressure gradient of 
Eq. (24) can be treated as a numerical (Neumann) boundary condition that allows the application 
of the divergence- free boundary condition u x + v y = 0. 

As discussed previously a finite-volume method is used for the discretization of the pressure 
equation. At a boundary point (i. j) this method is applied to a half cell, as shown in Fig. 2, 
surrounding the point. Let T a b c d be the boundary of the half cell. Along the surface boundary r a 6 



Figure 2: Half cell surrounding boundary point. 
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the pressure contribution to the integral of Eq. (18) can be expressed as 


pb pb 

/ Pxdy - Pydx = I Vp- 
J a J a 


n ds 



(25) 


where s = s(x, y ) is length along the boundary. The normal pressure gradient is given by Eq. (24) 
with the boundary condition u x + v y = 0 being imposed. On the boundary T sur f the normal 
pressure gradient becomes 


dp 

dn 


2 2 

■'EE Ti^dj dj V/c 

j = 1 k=l 

v\P'x(^ j xx “ 1 “ ^ yy ) “ 1 “ Tly(Vxx “ 1 “ V yy)\ 


(26) 


due to the no-slip condition. The integral of the normal pressure gradient is approximated as 
r b dp ( dp\ 

J ~dn^ S ^ ySny ^^ab = b / \fl x (Uxx T / U j yy')i,j T Tly(Vxx T Vyy'ji^^Sabi (27) 


where 


('U'xx T u yy)i,j 
i^xx T v yy)i,j 


\/2 E (■ u x Ay - u y Ax)r 

m= 1 
4 

^1/2 E (v x A y - VyAx) n 


m= 1 


with Ai /2 being the half cell area, and in the summations over the four edges of the half cell 
u x = — v y on r a fr. In the particular case of a flat plate (y = 0) the approximation of Eq. (27) 
becomes 


dp 

dn 


1,3 


dy. 


j A-S a b — v\(Vy)ij+ 1/2 ( \.Vy)i,j\As a t,A^i2 — l/ '{ v y)i,j+l/2 AXabA-^^- 


At the outflow boundary r on £ the x-momentum and y-momentum equations provide numerical 
boundary conditions for the calculation of the u and v velocity components. The boundary condition 
V • V = 0 is applied to the physical diffusion terms by substituting u x = — v y . The streamwise 
pressure gradient is approximated with second-order one-sided differencing. After updating the 
pressure along the boundary and solving the system 


( Qu 0 \ (6u\ = (R*\ 

V 0 Q h J \6v) \R h v ) ’ 


(28) 


the provisional updates u l+1 and v l+1 are computed. Then, the solution updates are determined 
by 

(;i;:)y::)y)r-w4 <»> 

where t x and t y are the Cartesian components of the unit tangent vector to the boundary T ou t. 

Suppose we define a sufficiently accurate value of boundary data as one consistent, to the order 
of approximation, with other data being specified. If a sufficiently accurate value for the tangential 



velocity (Vt) ou t at the boundary T ou t is not known, then a significant error in the solution at and 
near the boundary can be produced. That is, a numerical boundary layer can be created. Such an 
error can prevent attaining the desired order of accuracy simultaneously in the L\ and L 2 norms. 
This behavior was observed at the downstream boundary when computing parabola flows. 

An alternative to specifying ( Vt)out is to impose the condition that the streanrwise physical dif- 
fusion associated with the tangential momentum equation is zero. Then, the tangential momentum 
equation is solved for one velocity component, and the continuity equation is solved for the other 
velocity component. 

5 SOLUTION PROCEDURE 

The solution algorithm used for the steady, incompressible Navier-Stokes equations relies upon the 
factorizability of the corresponding discrete equations. In the present case factorizability necessarily 
implies a triangular discrete operator matrix. Such a matrix means that the elliptic operator A, 
which is associated with the pressure equation, and the advection-diffusion operator Q u can be 
inverted separately. Rapid inversion, in the sense of operation count, for each of the operators is 
achieved with a multigrid method. Line Gauss-Seidel relaxation is the smoother for the multigrid 
process. In the subsequent subsections the relaxation and multigrid schemes are described. 

5.1 Relaxation 

The relaxation scheme solves for u, v, and p along lines in each coordinate direction. This alternating 
line relaxation removes the influence of mesh cell aspect ratio on convergence rate. In order to 
describe the iteration matrix we consider for simplicity a rectangular domain with discrete points 
defined by {(xij, y t j) : 1 < i < m, 1 < j < n}. The iteration matrix for a relaxation sweep with 
line solves along lines of constant x is defined by 

B 2 C‘2 

^3 B : > C 3 

AA 4 A4 B4 C4 

AA m — 3 A m— 3 B m — 3 Cm — 3 

AA m - 2 A ni—2 B m — 2 C m — 2 

AA m — 1 

A m — 1 B m — 1 

where the elements of the change in the solution vector <5w and the residual function R are given 
by 


Swi = 

£ 

£ 

CO 

^Qi,n— 2 ^90, n— l) ) 

Svij 5p i:j ) T 

Ri = 

(g ,2 r *,3 • • • 

\ T 

^i,n— 2 G,n— 1 ) i I" i,j 

= (( r u)i,j ( r v)i,j i r p)i,j ) 


with i = 2, m — 1 and j = 2, n — 1. The submatrices AA, A, and C are (n — 2) x (n — 2) block 
diagonal matrices. The B submatrix is a block pentadiagonal matrix, where the diagonal blocks 
of the submatrices are 3x3 diagonal matrices. All of the other blocks of the submatrices are 
3x3 upper triangular matrices. The diagonal elements of the submatrices are the coefficients 
corresponding to the discrete approximations to the operators Q u and A. For convenience Eq. (30) 
can be rewritten as 

M8w = -R, (31) 
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where M is the iteration matrix and 

T T 

d'w = (i>w 2 <Jw 3 ••• <5w m _ 2 Sw m ,__ i ) , R=(R 2 R 3 ••• R m _2 R m -i) • 

By taking M to be the block lower triangular matrix resulting from ignoring the blocks above the 
principal diagonal (defined by the B matrices) , lexicographic line Gauss-Seidel is obtained with line 
solves corresponding to lines of constant x. The upper triangular form of the sub-blocks in each 
of the block matrices of M allows the solution of a tridiagonal system for the pressure on constant 
x lines followed by the solution of two pentadiagonal systems for the velocity components. Since 
the velocity on constant x lines is updated using a known pressure, the sequential updating of u 
and v is equivalent to collective relaxation. In a similar manner the iteration matrix for line solves 
along lines of constant y can be constructed. This basic part of the relaxation approach does not 
take into account the coupling between the equations that occurs because of boundary conditions 
and the stiffness of the equations in the vicinity of solid boundaries. As discussed subsequently the 
effects of coupling are treated with some additional relaxation at the boundaries. 

In numerous practical aerodynamic flows the cross flow physical diffusion dominates the stream- 
wise physical diffusion. If, for such flow problems, upwind differences are used for the advection 
operator of Eq. (3), then lexicographic line Gauss-Seidel relaxation is essentially equivalent to 
space marching of the advection terms. The advected error is effectively eliminated in one relax- 
ation sweep and the convergence rate becomes that of the Poisson equation for the pressure. This 
type of relaxation is sometimes called downstream relaxation. It is possible to get ideal multigrid 
convergence rates for the system because each component of the error is treated appropriately. 

5.2 Multigrid 

A sequence of grids G k , G k ~i , . . . , Go is used, where G*, is the finest and Go the coarsest. A 
straightforward full approximation storage (FAS) multigrid scheme is applied to the system of 
equations. In the following notation the indices (■ i,j ) are suppressed for convenience. Let L k be 
the discrete approximation to the operator matrix L and q^ be the solution on the fc-th grid. Also, 
let L k ~ i be the coarse-grid operator, l£' _1 be the fine-to-coarse grid restriction operator, and l£_ 1 
be the coarse-to-fine grid prolongation operator. If q*. is the current solution on grid k. the residual 
on this grid is r k = f k — L k q*,. This leads to the coarse-grid equation 

L k - iq fc _i = ffc_i = l£ _1 r fc + Z fc _i (/jT’qfc) . (32) 

After solving the coarse-grid equation for (\ k -\ , the fine-grid solution is corrected by 

q k <- qfc + i k -i (q*!-i - • (33) 

Equation (32) is solved by applying the same relaxation procedure that is used to solve the fine- 
grid equation. Multigrid is applied recursively to the coarse-grid equation. On the coarsest grid, 
many relaxation sweeps (typically 20 sweeps) are performed to ensure that the equation is solved 
completely. 

The intergrid transfer operators are standard ones (e.g., see [7] or [5]). A full weighting re- 
striction operator is used to transfer residuals. The transfer of the solution from a fine grid to a 
coarse grid is done with simple injection. Coarse-grid corrections are transferred with a bilinear 
interpolation operator, which is the adjoint of the restriction operator. A conventional Recycle is 
used to execute the multigrid process. 
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6 SKIN FRICTION 


For wall-bounded shear layers accurate prediction of surface skin friction is generally considered 
an essential requirement of an algorithm for computing viscous flows. In the case of Cartesian 
coordinates the skin friction is defined by 


C f = 


“2t w 

~Re' 


T„, = U. 


y, 


(34) 


where the velocity component parallel to the surface u is nondimensionalized by the free-stream 
velocity and the coordinate normal to the surface y is scaled by the reference length for the Reynolds 
number Re. Thus, the surface shear stress t w is determined by the normal gradient of the velocity 
parallel to the solid surface. In the case of curvilinear coordinates (£, rj) the wall shear stress can 
be computed from 

Vr/ 

Tw = W t • n, n = — — , 

I Vr/| 

where Vt is the velcity component parallel to the surface and n is the unit normal to the surface. 
Assume that a constant y curve coincides with the solid surface. Also, let J be the Jacobian for 
the transformation (x,y) — ► (£, 77 ). The inner product of Eq. (35) gives 



Tw 


{x\ + yf) 1/2 

F 1 


(VtF 


(36) 


and J 1 = x^yn - x v y^. 

In cell-centered finite-volume schemes the wall shear stress is often computed with a two-point 
one-sided difference for V) divided the normal distance from the first interior point to the wall. A 
Taylor series expansion shows that this approximation is second order only when the velocity profile 
is linear near the surface. With the present incompressible scheme, which uses finite-differencing 
for the advection terms, such an approximation was always first order. Even in the case of Hiemenz 
flow, which is planar stagnation flow without surface curvature, the behavior of this approximation 
was first order. In the numerical computation of t w the quantity (Vt) v is discretized with a second- 
order one-sided difference. The transformation derivatives x v and y v are also approximated with 
one-sided differences. The derivatives x % and are approximated with central differencing. 

The skin-friction drag coefficient is determined by 



Cfds , 


(37) 


where s is distance along the surface boundary T sur f. The trapezoidal rule is used for the numerical 
evaluation of the integral in Eq. (37). With the second-order approximation for (V)), ; the computed 
(Cd) f for Hiemenz and circular-arc stagnation flows, as well as parabola flow, was found to be 
second order with mesh refinement. 


7 NUMERICAL RESULTS 

The numerical scheme described in the previous sections was applied to three incompressible, 
viscous flow problems. As an initial evaluation of the scheme, high Reynolds number flow past a 
flat plate at zero incidence was considered. Then the effect of geometric curvature and Reynolds 
number on the behavior of the present scheme was examined by considering flow around a parabola 
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x/L 


Figure 3: Domain of flat plate flow with 33 x 33 grid. 

(representing an airfoil leading edge). For the final problem the scheme was applied to flow over 
a symmetric Karman-Trefftz airfoil. Computational results are given for two different angles of 
attack. 

7.1 Flat Plate 

The computational domain for the flat plate simulation is displayed in Fig. 3. The plate begins at 
x/L = 0, the inflow boundary is at x/L = 0.2, and the outer boundary is at x/L = 1. The upper 
boundary is located one plate length away from the plate (y/L = 1). At the inflow boundary, 
the Blasius (self-similar) solution was prescribed. For the upper and downstream boundaries the 
pressure was set to the free-stream value p^, and the velocity components tangential to these 
outflow boundaries were defined according to the Blasius solution. The no-slip condition was 
imposed on the plate (x/L = 0.2 to x/L = 1). 

The finest grid used for the flat plate calculations consisted of512x512 cells. The grid spacing 
in the x direction was uniform. In order to resolve the boundary layer on the plate the grid was 
clustered at y/L = 0 and stretched geometrically to y/L = 1 (see Fig. 3). For this grid the 
minimum spacing in the y direction was 0.000375, and the stretching factor was 1.0053. A series of 
nested coarse grids was obtained by coarsening the fine grids by a factor of two in each coordinate 
direction. In all cases shown below, the coarsest grid was 16 x 16 cells. 

Line Gauss-Seidel relaxation was used for the computations. The complete relaxation process 
involved one sweep with vertical line solves and one sweep with horizontal line solves, along with 
some additional work near the plate. The relatively small amount of additional work is needed 
due to the coupling of the flow equations near the plate. Moreover, the points for j < 3 were 
relaxed with 3-5 sweeps of horizontal solves. For the vertical line solves relaxation was initiated 
at the inflow boundary (x/L = 0.2) and proceeded to the outflow boundary (x/L = 1). With 
vertical solves underrelaxation was applied at 5-10 points near the plate. The underrelaxation 
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factor was (1 + /3) -1 , where (3 was between 0.5 and 1.0. In the case of the horizontal line solves 
relaxation started at the outer boundary ( y/L = 1) and continued to the inner boundary (y/L = 0). 
A W(l, 1) multigrid cycle was used; that is, one complete relaxation process was performed on each 
grid before restricting to the coarse grid, and one complete relaxation process was performed after 
the coarse- grid correction was added to the fine-grid solution. 

The Reynolds number of the flow past the flat plate was 10, 000. Figures 4 and 5 show the 
variation of the tangential velocity component u and the transverse velocity component v with the 
scaled normal coordinate r/. Both velocity components are nondimensionalized by the boundary- 
layer edge velocity (u e ). The scaled v is multiplied by ReJ , where Re x is the Reynolds number 
based on x, to obtain the self-similar solution, and the coordinate q = 5 y/8, where 6 is the thickness 
of the boundary layer. The computed velocity profiles are at the midplate location. Starting with 
the 32 x 32 cell grid, which has only seven points in the boundary layer, there is excellent agreement 
with the classical Blasius tangential velocity. There is good agreement with the Blasius transverse 
velocity with just seven points in the boundary layer. 

In Figs. 6 and 7 the convergence behavior of the scheme for the flat plate flow is shown. The 
L -2 norm of the residuals for the pressure and momentum equations is given for each W(l, 1) cycle 
in Fig. 6. In the full multigrid (FMG) process a single grid calculation was performed on the 
17 x 17 grid. The multigrid algorithm was applied at each successive level of the FMG process. 
For most of the grids the residuals are reduced more than 6 orders of magnitude in 10 cycles. The 
convergence rates on the finest grid (513 x 513) for the three flow equations are 0.12 to 0.15. On 
the two finest grids the residuals are essentially reduced to machine zero. According to local mode 
analysis (Brandt[13]) the smoothing factor for relaxation with line solves is 5 -1 / 2 . Assuming only 
radial line solves in the W(l,l) cycle, the residual reduction per cycle would be 0.2. Due to the 
azimuthal line solves a faster convergence rate is attained. 

Figure 7 shows the variation of the ratio of algebraic error to discretization error for the skin- 
friction drag coefficient of the plate. The skin-friction Cf is calculated from x/L = 0 to x/L = 0.2 
by assuming it has an inverse square root behavior in distance from the leading edge of the plate 
(see Thomas, Diskin, and Brandt [9]). Numerical integration with the trapezoidal rule of Cf over 
the entire plate yields the drag coefficient. For all levels of grid refinement the algebraic error is 
reduced below the truncation error in a single multigrid cycle. Such performance of a scheme may 
be interpreted as an alternative definition of an optimally convergent scheme. 

As indicated in Figs. 8 and 9, the solution on the finest grid can be obtained in a FMG process 
involving only one multigrid cycle on each grid considered in Fig. 6. The behavior of surface skin 
friction with mesh refinement is displayed in Fig. 10. In this figure N y denotes the number of 
points in the y-direction. Also, the zero subscript refers to a mesh with 97 points in the normal 
direction and a stretching factor of 1.03 when the upper boundary of the domain is at y m ax/L = 1.0. 
Calculated results are given for two different locations of the upper boundary. The linear variation 
of the skin-friction ratio C / / (C /) Blasius with the mesh spacing squared demonstrates second-order 
accuracy for the skin friction. 
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Figure 4: Tangential velocity profiles at midplate location for laminar flow over flat plate 
(. Re L = 10,000). 



(v/u e ) * sqrt(Re x ) 


Figure 5: Transverse velocity profiles at midplate location for laminar flow over flat plate 
(. Re L = 10,000). 
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17x17 33x33 65 x65 129 x129 257x257 513x513 



Figure 6: Convergence behavior for laminar flow over flat plate (ReL = 10,000). Residual 
symbols: Circle for pressure, square for x-niomentum, triangle for ^/-momentum. 



Figure 7: Convergence behavior of drag coefficient for laminar flow over flat plate ( Re^ = 
10 , 000 ). 
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Figure 8: Variation with multigrid cycles of tangential velocity profile at midplate location 
for laminar flow over flat plate (Re^ = 10, 000, 129 x 129 grid). 



Figure 9: Variation with multigrid cycles of transverse velocity profile at midplate location 
for laminar flow over flat plate (Rei = 10, 000, 129 x 129 grid). 
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Figure 10: Behavior with mesh refinement of surface skin friction for laminar flow over flat 
plate ( ReL = 10,000). 
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x/R 

Figure 11: Coarse grid (33 x 17) for laminar flow over parabola. 

7.2 Parabola 

The flow over a parabola was determined with a conformal mapping. For the parabola, the leading 
edge radius of curvature R is the reference length. In the mapping a rectangular region in the 
complex z-plane was mapped to a region in a complex ((-plane with inner (solid surface) and outer 
(inflow) boundaries defined by parabolas. The surface length of the inner boundary is 2, and the 
distance between the vertices of the inner and outer parabolas is 1. The physical problem in the z- 
plane is planar stagnation (Hiemenz) flow. In the computations the velocity at the outer boundary 
is specified according to the complex velocity in the ((-plane. At the outflow boundary the pressure 
and the velocity tangent to the boundary are specified, resulting in a well-posed problem for the 
full Navier-Stokes equations. The pressure is the free-stream value, and the tangential velocity is 
obtained from the complex velocity in the ((-plane. 

Laminar flow over the parabola at three different Reynolds numbers was considered. The 
Reynolds numbers based on radius of curvature were ReR = 25, Rcr = 100, and Ren = 312.5, 
which correspond to Reynolds numbers based on unit airfoil chord of Re c = 800, Re c = 3200, and 
Re c = 10000, respectively. In Fig. 11a coarse grid for the calculations is shown. Figure 12 presents 
the variation in the error of the continuity equation with mesh refinement. Second-order accuracy 
of the scheme is essentially realized on the 129 x 65 grid for all Reynolds numbers. In Fig. 13 the 
change in surface skin-friction error, at the location midway between the vertex of the parabola 
and the upper outflow boundary, with mesh density is shown. Second-order accuracy begins on the 
65 x 33 grid. Numerical experiments indicate that the principal parameter affecting the onset of 
second-order behavior is the curvature at the leading edge of the parabola. Moreover, as the radius 
of curvature is increased, the mesh density for the onset of design accuracy decreases. 

In Table 1 second-order accuracy is also demonstrated for the drag coefficient ( Cd ) of the 
parabola at each ReR. The value of the Cd is given for both one multigrid cycle (FMG-1) and ten 


18 





Figure 12: Behavior of continuity with mesh refinement and Reynolds number for laminar 
flow over parabola. Mesh range: 17 x 9 — > 513 x 257. 



Figure 13: Variation of skin friction error with mesh refinement and Reynolds number for 
laminar flow over parabola. Mesh range: 17 x 9 — > 513 x 257, Axial location: x/R = 16.59. 
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Ren 

N x X Ny 

C d (FMG-1) 

C d (FMG-10) 

|e a |/|e d | : C d 
(FMG-1) 

25 

17 x 9 

1.070609 

1.030401 

0.029153 


33 x 17 

1.489693 

1.535983 

0.283019 


65 x 33 

1.603028 

1.617511 

0.288438 


129 x 65 

1.641015 

1.643881 

0.234462 


257 x 129 

1.649742 

1.651158 

0.404654 


513 x 257 

1.649895 

1.652720 

0.844388 

100 

17 x 9 

0.466985 

0.445406 

0.041222 


33 x 17 

0.552725 

0.677579 

0.512038 


65 x 33 

0.725510 

0.758993 

0.469649 


129 x 65 

0.783574 

0.787396 

0.288910 


257 x 129 

0.794248 

0.794577 

0.128796 


513 x 257 

0.795718 

0.796246 

0.487271 

312.5 

17 x 9 

0.224233 

0.208086 

0.059588 


33 x 17 

0.230412 

0.344023 

0.534571 


65 x 33 

0.363082 

0.397327 

0.427955 


129 x 65 

0.423235 

0.431337 

0.407791 


257 x 129 

0.440050 

0.440120 

0.022917 


513 x 257 

0.441890 

0.442357 

0.384757 


Table 1: Drag coefficient at different Reynolds numbers for laminar flow over parabola. 
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Figure 14: Effect of Reynolds number on solution convergence for laminar flow over parabola. 
Residual symbols: circle for Ren = 25, square for Ren = 100, triangle for Ren = 312.5. 
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cycles (FMG-10). In addition, the ratio of the algebraic error to the discretization error, |e a |/|e<i|, 
is presented for each grid. The algebraic error is reduced below the truncation error in one cycle 
for all grids and all Reynolds numbers. Thus, for this flow over a body with curvature, the present 
scheme satisfies the alternative definition of an optimally convergent scheme. 

In Fig. 14 the effect of Reynolds number on the convergence behavior for the parabola flow is 
shown. The residual is for the pressure equation. Convergence behavior is similar for the higher 
mesh densities (i.e., 129 x 65 and greater) at each Reynolds number, and rapid convergence with 
at least a rate of 0.2 per cycle is attained for these mesh sizes. As the mesh density decreases (i.e., 
65 x 33 and smaller) , the convergence rate becomes slower and the influence of Reynolds number on 
the rate becomes greater. Such coarse grid slowdown has been observed for inviscid flows (Roberts 
and Swanson [11]) as well as viscous flows being considered here. One possible explanation for this 
behavior is the larger discretization error on coarser grids associated with the relaxation process. 
Thus, whenever the resolution on the coarser grids in the FMG procedure is not adequate, one 
would anticipate convergence slowdown. 

7.3 Airfoil 

Solutions were obtained for incompressible, viscous flow around a symmetric Karman-Trefftz airfoil 
at zero and two degrees angle of attack. The Reynolds number based on chord ( Re c ) was varied 
from 200 to 800. A Karman-Trefftz airfoil was generated from a cylinder by a conformal mapping 
[19]. The trailing-edge angle is 10°, and the thickness ratio is approximately 15 percent. The airfoil 
flow was solved on a finite domain. At inflow points along the outer boundary the free-stream 
velocity components were specified. For outflow points the free-stream pressure was prescribed, 
and the streamwise diffusion for the tangential momentum equation was taken to be zero. 

The present scheme is compared with Runge-Kutta (R-K) and implicit upwind based multigrid 
schemes. These other two schemes solve the compressible Navier-Stokes equations, and they do 
not include preconditioning. Thus, as the Mach number goes to zero, these schemes can experience 
inaccurate discretization and slow convergence, or even divergence. In order to have representative 
convergence of the schemes, and at the same time have the Mach number sufficiently low for 
comparison purposes, the free-stream Mach number (M^) for the calculations with the R-K and 
implicit upwind schemes was taken to be 0.1. The downstream outflow boundary pressure p e was 
specified in the calculations with the compressible flow solvers. Thus, to simulate the incompressible 
flow and be able to compare with the solution from the incompressible flow solver, it was necessary 
to increment the pressure p e . For example, with the outer boundary at three chords, the modified 
pressure p* e = p e + A p = 1.00015 when M ^ = 0.1, and A p = 0.00015. As — > 0, then A p — > 0. 

To facilitate the resolution of the near and far wake regions, especially in the case of laminar 
flow where the mixing rate due to the physical diffusion is slow, a C-type mesh was generated for 
the airfoil calculations. A hyperbolic grid generator was used to create the mesh. The near field 
part of the grid is depicted in Fig. 15. The outer boundary of the domain is roughly 3 chords 
from the airfoil. The coarsest grid in the grid sequence for the multigrid solver contained 16 x 8 
cells. On each grid a relaxation sweep started at the radial line emanating from the airfoil leading 
edge, proceeded over the upper half of the domain, and then over the lower half of the domain. 
The relaxation sweep with radial solves was followed by a relaxation sweep with azimuthal solves 
that started at the outer C-type boundary. Then, 3-5 additional sweeps with azimuthal solves were 
performed near the solid surface boundary (j < 3). This is done to reduce the residuals near the 
boundary so that they are comparable in size to residuals in the interior. 

In Fig. 16 the convergence histories for this nonlifting case are presented. The average rate 
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Figure 15: Near field of 128 x 64 cell C-type grid for computing laminar flow over Karman- 
Trefftz airfoil (. Re c = 200). 


17x9 33 x17 65 x33 129x65 257x129 513x257 



Figure 16: Convergence behavior for laminar flow over nonlifting Karman-Trefftz airfoil 
(' a = 0°, Re c = 200). Residual symbols: Circle for pressure, square for ^-momentum, 
triangle for {/-momentum. 
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Scheme 

MG Cycle 

WU/cycle 

Tridiag Solves/ 
N X Ny 

N cycle 

\\Res\\ 2 < 10- 4 

INCOMP 

W(l,l) 

10.7 

20 

5 

5 stage R-K 

W(1,0) 

12 

36 

124 

Implicit upwind 

W(1,0) 

4 

12 

485 


Table 2: Computational work for incompressible, 5-stage R-K, and implicit upwind schemes. 

of reduction of the residual on the 513 x 257 grid is about 0.11 for the pressure and momentum 
equations. There is again somewhat slower convergence on the coarser meshes, with about five 
orders of magnitude reduction of the residuals on the 65 x 33 grid. 

In order to estimate the relative work required by the present, R-K, and implicit upwind 
schemes, we consider an operation count based on residual evaluations, tridiagonal solves, and 
the number of cycles to achieve a specific reduction in residuals. First we define the work as- 
sociated with residual evaluations. For a W (m,n) multigrid cycle, there are m + n + 1 residual 
evaluations per cycle on the fine grid. Let WU represent work units. A work unit is often defined 
to be one relaxation sweep of the finest grid. Thus, a WU can be considered as the cost of a vector 
residual evaluation on the finest grid. In this assessment we neglect the cost of interpolating the 
residuals and solutions between the grid levels. For a W (m,n) cycle, we have that 

wu r, , x ,w n 1 1 

H^"[' ; < m+ ") + l](l + 2 + 4 + "') 

= 2[k(m + n) + 1], (38) 

where k denotes the number of sweeps on each grid that require a vector residual evaluation. If 
line solves are considered in each coordinate direction, then k = 2. In the case of a five stage R-K 
scheme, k = 5. The WU per cycle for the current scheme with a W(l,l) cycle that includes one 
radial relaxation sweep and one azimuthal relaxation sweep (neglecting additional work at surface 
boundary with azimuthal solves) is 10. 

When a finite-volume formulation is used for the discretization of the flow equations and al- 
ternating line Gauss-Seidel relaxation is applied as in the scheme of Swanson [12], an additional 
residual evaluation is incurred due to the Gauss-Seidel updating of the solution. Since a finite- 
difference discretization is used for the momentum equations in the present scheme, there is only a 
1/3 additional vector residual calculation (due to the pressure equation) for each relaxation sweep. 
If we include that in the operation count, the WU per cycle for the present scheme would be in- 
creased to approximately 10.7. Both the R-K and implicit upwind schemes apply a W(1,0) cycle, 
and the WU per cycle for these schemes is given in Table 2. 

For the three schemes there is a difference in computational effort resulting from scalar tridi- 
agonal and pentadiagonal inversions. As indicated previously each relaxation sweep of the present 
scheme requires one tridiagonal solve and two pentadiagonal solves (corresponding to the three flow 
equations) for the grid lines not in the direction of the relaxation. Since the operation count for a 
pentadiagonal inversion is about two times that for a tridiagonal inversion, the work required on 
the fine grid with m + n = 2 and k = 2 is roughly proportional to that for 20 tri diagonal solves. 
Since the present scheme is solving three equations, we assume, for an equitable comparison of 
operation count, that all schemes are solving three equations. For the R-K scheme there is a scalar 
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Scheme 

res / (N 0 p S ) res 

(N 0 ps) tri 

Work 

INCOMP 

54 

200 

8296 

5 stage R-K 

1488 

8928 

255936 

Implicit upwind 

1940 

11640 

333680 


Table 3: Comparison of work required by three Navier-Stokes schemes (grid: 513 x 257). 

tridiagonal inversion for each flow equation, each coordinate direction, and each of the five stages 
due to the residual smoothing procedure used to extend stability. In addition, there is similar work 
performed on the prolongated coarse grid corrections. Thus, the work due to this part of the scheme 
is proportional to 36 tridiagonal solves. Thus, the R-K scheme considered here [20] requires almost 
two times as many tridiagonal inversions as the present scheme. It should also be pointed out that 
the R-K scheme used here is a (5,3) scheme, which means that the complete residual, including 
viscous and numerical dissipation terms, is evaluated only on three stages. This factor has not been 
taken into account. With the implicit upwind scheme there are six scalar tridiagonal inversions 
because a diagonalized version of the scheme is used. In the diagonalization procedure there are 
three matrix- vector multiplies, and for a system of three equations, these correspond to roughly six 
tridiagonal solves. Thus, the implicit upwind scheme requires about 12 scalar tri diagonal solves. 

We now sum the work contributions due to residual evaluations and tridiagonal solves. Let W res 
and Wt r i denote the work associated with residual evaluations and tri diagonal solves, respectively. 
Then, the total work per grid point required by a scheme can be approximated by 

Work = W res + Wtri (39) 

where 

W res = (WU/cycle) ( N ops ) res N cydes , (40) 

Wtri = 2 N tr i (N ops ) tri N cycles i (41 ) 

with N ops representing the number of operations needed to evaluate the residual vector function, 
Ntri being the number of tridiagonal solves per total number of points (N x N y ), ( N ops )t r i denoting 
the number of operations in a scalar tridiagonal solve (3 adds, 3 multiplies, and 2 divides for a 
total of 8 operations), and N cyc ies being the number of cycles needed to reduce the residuals a 
predetermined number of decades. The factor of two in Eq. (41) is due to the multigrid process. 
According to [22] the residual (based on central differencing with numerical dissipation) of the 
diagonalized implicit algorithm requires 124 operations (i.e., 48 adds, 72 multiplies, 4 divides). 
Then, an estimate for the ratio of ( N ops ) res to ( N ops ) tr i is 16. If we assume that ( N ops ) res is about 
the same for all schemes, the contributions of both W res and Wtri to the Work (see Table 3) 
indicate that the present scheme is more than 30 times faster than the other schemes for the finest 
mesh (i.e., 513 x 257). For a mesh density of 129 x 65 the speedup of the present scheme exceeds an 
order of magnitude. This improvement in the computational effort relative to the R-K and implicit 
upwind schemes is also evident in the variation of the drag coefficient [Cy) with multigrid cycles 
(see [12]). 

The total drag coefficient (Cy), along with the pressure and skin-friction contributions (( Cy) p 
and ( Cy ) f, respectively) computed with the three schemes on the 513 x 257 grid are given in Table 4. 
The maximum variation in Cy divided by the INCOMP scheme value is 2.2 percent. There is an 
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Scheme 

(C d ) p 

(C d )f 

c d 

INCOMP 

0.074780 

0.222239 

0.297019 

5 stage R-K 

0.076662 

0.226987 

0.303649 

Implicit upwind 

0.071351 

0.227320 

0.298670 


Table 4: Drag coefficients for laminar flow over Karman-Trefftz airfoil computed with three 
different schemes (a = 0°, Re c = 200, grid: 513 x 257). 

unexpected difference between the ( Cd) p for the implicit upwind scheme and the other schemes. 
The error behavior of the drag coefficients with mesh refinement was determined. For all algorithms 
the variation in ( C d )p error (based on a Richardson extrapolated value) with mesh refinement is at 
least second order. The error in (Cd) f obtained with the R-K and implicit upwind schemes does not 
decrease in a monotonic manner with reduction in mesh spacing. Although there is a monotonic 
decrease in ( C d ) / error for the INCOMP scheme, the error is only first order. Thus, none of the 
algorithms exhibit second-order accuracy for the total drag error. 

In Figs. 17-22 the pressure and skin-friction distributions computed with the present scheme on 
various grid densities are compared with the fine grid results obtained with the R-K scheme (where 
the free-stream Mach number was set to 0.1). There is generally very good agreement starting with 
the 129 x 65 grid for both the nonlifting and lifting flows. Figures 23 and 24 show the variation 
in the error of pressure and skin-friction near the location of peak skin friction. The error of each 
quantity is determined from the corresponding Richardson extrapolation value. Both the pressure 
and skin friction exhibit, in this region of strong flow gradients, essentially second-order behavior. 

The behavior of solution error with respect to mesh refinement is displayed in Fig. 25 and 
Fig. 26. The solution on the 513 x 257 grid is used to determine the error. Figure 25 shows that 
the L -2 norms of the errors, excluding boundary points, in the solution variables u, v, and p exhibit 
second-order accuracy. In addition, this figure reveals that the error in the continuity equation 
decreases in a first-order manner. The derivatives in the continuity equation are approximated 
with central differencing. Error variations in the flow quantities along the wake line ( x/c > 1, 
y/c = 0) are presented in Fig. 26. The pressure error diminishes with second-order accuracy. As 
the streamwise grid stretching factor (equal to 1.075 on the finest mesh) in the wake decreases with 
mesh refinement, the error in the velocity component u is approaching second-order behavior. In 
the case of the continuity equation, the error appears to be going to a constant as the mesh is 
refined. A mass conservation error is generated at the airfoil trailing edge, and this error is being 
transported by advection and diffusion as the flow proceeds downstream. The occurrence of such 
an error certainly indicates that an improvement is required in the accuracy of the discretization of 
the pressure equation (i.e., boundary condition) at the trailing edge. At the same time, it should 
be emphasized that the solution on the airfoil does not appear to be significantly compromised, as 
comparisons with solutions computed with other methods have demonstrated. 

Pressure and velocity contours for the lifting solution on the 513 x 257 C-type grid are depicted 
in Figs. 27 and 28. The thickness of the attached boundary layer is about 0.25 airfoil chord, and 
there are about 70 points in the boundary layer at the midchord location. Slow diffusion of the 
airfoil wake is evident. 

As indicated previously the outer boundary (OB) for the initial set of calculations was located 
only about 3 chords (approximately 12 boundary-layer thicknesses) away from the surface of the 
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Re c 

(C d ) P 

{Cd)f 

c d 

200 

0.070366 

0.211714 

0.282080 

400 

0.056302 

0.141830 

0.198132 

600 

0.049628 

0.112108 

0.161736 

800 

0.045483 

0.094836 

0.140319 


Table 5: Variation of drag coefficient with Reynolds number for laminar flow over airfoil 
(Karman-Trefftz, a = 0°, grid: 641 x 449). 


Re c 

(C d ) P 

{Cd)f 

c d 

200 

0.070210 

0.211399 

0.281609 

400 

0.056185 

0.141690 

0.197875 

600 

0.049527 

0.112031 

0.161558 

800 

0.045390 

0.094794 

0.140184 


Table 6: Variation of drag coefficient with Reynolds number for laminar flow over airfoil 
(Karman-Trefftz, a = 0°, grid: 897 x 513). 

airfoil. Thus, a series of nonlifting computations were performed to assess the influence of the OB 
location on the solutions. The OB location was varied between 3 chords and 20 chords. For each 
different location of the boundary the same grid stretching factor was maintained in the directions 
normal to the airfoil surface and downstream of the airfoil trailing edge. The effect of the OB 
distance, which is designated by R, on the skin-friction and total drag coefficients is displayed in 
Fig. 29. The error in the two drag coefficients is based on values determined from Richardson 
extrapolation. There is only a 1.6 percent error in total drag when the OB is located at three 
chords. 

The performance of the present scheme for the airfoil flow is now considered for different 
Reynolds numbers (Re c ). In Table 5 the variation of the pressure, skin-friction, and total drag 
coefficients with Re c are given. These nonlifting solutions were computed with a mesh density of 
641 x 449 and an OB location of 20 chords. This set of calculations was repeated with the tangen- 
tial mesh spacing on the airfoil decreased by a factor of two and a reduction in the normal mesh 
spacing, resulting in a mesh density of 897 x 513. The drag coefficients for these computations are 
presented in Table 6. Differences between these two sets of results appear in the fourth decimal 
place. In Fig. 30 the convergence behavior on the 641 x 449 grid for Re c = 200 is shown. The rates 
of convergence on the sequence of grids are essentially the same as those obtained for the sequence 
associated with the 513 x 257 grid having the OB at the 3 chord location. Convergence histories 
for Re c = 600 and mesh densities of 641 x 449 and 897 x 513 are displayed in Figs. 31 and 32. The 
convergence rates on both mesh sequences are similar. On the 641 x 449 grid the convergence rate 
(average residual reduction rate) is about 0.13 per cycle. For the 897 x 513 grid the convergence 
rate is slightly slower (0.15 per cycle). 

Figure 33 shows the logarithmic variation of the ratio of algebraic error to discretization error 
for the skin-friction drag coefficient of the airfoil. The present scheme nearly exhibits optimal 
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Re c 

a 

(C d ) P 

(C d )f 

c d 

200 

0.089128 

0.072462 

0.210588 

0.283050 

400 

0.084476 

0.058133 

0.140861 

0.198994 

600 

0.079284 

0.051275 

0.111242 

0.162517 


Table 7: Variation of drag coefficient with Reynolds number for laminar flow over lifting 
airfoil (Karman-Trefftz, a = 2°, grid: 641 x 449). 

multigrid behavior, according to the alternative definition given previously. Only on the finest grid 
is the algebraic error not driven below the truncation error of the scheme in one multigrid cycle. 

Solutions for the lifting Karman-Trefftz airfoil were also obtained on the 641 x 449 grid with the 
OB at 20 chords. Computed drag coefficients corresponding to Re c values of 200, 400, and 600 are 
given in Table 7. Convergence difficulties were encountered at Re c = 800. A subsequent calculation 
done with the R-K scheme indicated flow separation on the airfoil upper surface near the trailing 
edge at this Reynolds number. The issues related to separated flows have not been considered in 
this report. Some of these issues, such as poor coarse-grid approximation, have been addressed 
by Brandt and Yavneh [8]. In addition, fast convergence has been obtained for the driven cavity 
problem, provided the Reynolds number of the flow is sufficiently low (i.e., Re < 1000). However, 
further effort is required to realize a second-order accurate scheme with optimal convergence when 
solving steady flow problems over a large Re range that have streamwise recirculation regions. 

An effort to consider the effect on convergence rate of refinement in normal airfoil surface grid 
spacing and grid stretching was made. Reduction in the normal spacing required an increase in 
the magnitude of the underrelaxation parameter /3, which is used near the airfoil surface when 
relaxing with radial lines. A doubling of the standard value was necessary to achieve stability 
and convergence. This increase in (3 introduced an unsatisfactory mesh sensitivity with respect 
to convergence. In addition, larger (3 values result in slower convergence. A possible remedy for 
this difficulty is to solve the flow equations fully coupled near the surface. Then there will be no 
requirement for an underrelaxation factor, and the robustness of the scheme should be significantly 
enhanced. 


8 CONCLUDING REMARKS 

A multigrid algorithm with essentially O(N) behavior has been developed for the steady, incom- 
pressible Navier-Stokes equations. The momentum and pressure equations have been discretized 
with conventional finite- difference and finite-volume methods, respectively. A discrete approxima- 
tion to the normal component of the momentum equation has been used as a vehicle for applying 
the divergence-free boundary condition for the pressure equation. Solutions have been obtained 
for laminar flow past a flat plate with essentially textbook multigrid efficiency. Furthermore, the 
multigrid method has been used to solve viscous flows with curvature effects. These geometries 
include a parabola and a Karman-Trefftz airfoil. Fast convergence has been obtained for each of 
these flows. With the present scheme more than an order of magnitude reduction in computational 
effort has been achieved when compared with R-K and implicit upwind based multigrid schemes. 

Solutions have been calculated for both nonlifting and lifting airfoil flows with Reynolds number 
(Re c ) varying between 200 and 800. Tables for aerodynamic coefficients have been presented. 




For the lifting case with Re c = 800 separation occurred and convergence was not possible. The 
capability to compute separated flows was outside the scope of this report. 

It has been demonstrated that the dependent variables of the airfoil calculations exhibit a global 
second-order accuracy. Furthermore, it has been shown that the dependent variables have an order 
property in the wake region. However, these computations also indicate that an error is created in 
mass conservation at the airfoil trailing edge that is convected and diffused downstream. This error 
compromises the order property of the continuity equation in the wake region. Although there 
appears to be no significant effect on the solution in the vicinity of the airfoil, such behavior in the 
continuity equation certainly is not considered acceptable. Preliminary numerical results indicate 
that a stronger enforcement of continuity near the trailing edge can provide the necessary remedy. 

The underrelaxation factor (/5) , which is used in conjunction with radial line solves at grid 
points near a solid surface boundary, can cause a slowdown in convergence. Such behavior has 
occurred when the normal mesh spacing at the surface becomes sufficiently small that the coupling 
of the flow equations requires a significant increase (i.e., a factor of two) in the standard (3 to 
maintain stability. This difficulty can be removed by solving the flow equations fully coupled near 
the surface. Then there would be no requirement for underrelaxation. 

While fast convergence has been obtained on fine grids, the convergence on relatively coarse 
meshes in the FMG process has been somewhat slow. The sources of slower convergence rates on 
coarse grids need to be completely identified. One possible factor that can affect these convergence 
rates is the validity of the principal linearization on the coarse grids. That is, terms considered 
subprincipal for the relaxation process may have become principal. Another possible factor affecting 
these rates is the accuracy of the numerical boundary conditions as well as the choice of boundary 
conditions. 

Further studies of the present scheme are also required to clearly delineate the effects of mesh 
aspect ratio and grid stretching on convergence rate, and thus, on the discrete factorizability of the 
algorithm. While significant gains in computational efficiency have been achieved for laminar flow, 
much greater gains (i.e., two orders of magnitude) are anticipated with the extension of factorizable 
schemes to allow computation of turbulent flows. 

One viewpoint of the present scheme is that it represents a building block in a general class of 
algorithms for solving complex flows. For separated and turbulent flows the INCOMP scheme will 
have to perform in concert with other building blocks in order to attain optimal multigrid efficiency. 
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Surface pressure distribution for laminar flow over nonlifting Karman-Trefftz 
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Figure 18: Surface skin-friction distribution for laminar flow over nonlifting Karman-Trefftz 
airfoil (a = 0°, Re c = 200). 






Figure 19: Blowup of Surface skin-friction distribution for laminar over nonlifting Karman- 
Trefftz airfoil (a = 0°, Re c = 200). 



Figure 20: Surface pressure distribution for laminar flow over lifting Karman-Trefftz airfoil 
(a = 2°, Re c = 200). 
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Figure 21: Surface skin- friction distribution for laminar flow over lifting Karman-Trefftz 
airfoil (a = 2°, Re c = 200). 



Figure 22: Blownp of Surface skin-friction distribution for laminar over lifting Karman- 
Trefftz airfoil (a = 2°, Re c = 200). 
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Figure 23: Variation in pressure coefficient error with mesh refinement for laminar flow over 
Karman-Trefftz airfoil. Mesh sequence: 33 x 17 — > 513 x 257; x/c = 0.021; a = 0°,2°; 
Re c = 200. 



Figure 24: Variation in skin-friction coefficient error with mesh refinement for laminar flow 
over Karman-Trefftz airfoil. Mesh sequence: 33 x 17 — » 513 x 257; x/c = 0.021; a = 0°,2°; 
Re c = 200. 
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Figure 25: Behavior of error in u, v , p, and continuity with mesh refinement for laminar flow 
over Karman-Trefftz airfoil. Mesh sequence: 17 x 9 — » 257 x 129; a = 0°; Re c = 200. 
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Figure 26: Wake line behavior of L 2 error in u, p, and continuity with mesh refinement for 
laminar flow over Karman-Trefftz airfoil. Mesh sequence: 17 x 9 — > 257 x 129; a = 0°; 
Re c = 200. 
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Figure 27: Pressure contours for laminar flow over lifting Karman-Trefftz airfoil (a = 2°, 
Re c = 200). 



Figure 28: Velocity contours for laminar flow over lifting Karman-Trefftz airfoil (a = 2°, 
Re c = 200). 
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Figure 29: Error in skin-friction and total drag coefficient due to outer boundary location 
(Karman-Trefftz airfoil, a = 0°, Re c = 200). 


21 X 15 41 X 29 81 x57 161 x113 321 X 225 641 X 449 



Figure 30: Convergence behavior for laminar flow over nonlifting Karman-Trefftz airfoil 
(' a = 0°, Re c = 200, OB = 20 chords). Residual symbols: Circle for pressure, square for 
x-momentum, triangle for ^/-momentum. 
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21 X 15 41 X 29 81 x57 161 x113 321 X 225 641 X 449 



Figure 31: Convergence behavior for laminar flow over nonlifting Karman-Trefftz airfoil 
(i a = 0°, Re c = 600, OB = 20 chords). Residual symbols: Circle for pressure, square for 
x-momentum, triangle for ^/-momentum. 


29 x17 57x33 113x65 225x129 449x257 897x513 



Figure 32: Convergence behavior for laminar flow over nonlifting Karman-Trefftz airfoil 
(a = 0°, Re c = 600, OB = 20 chords). Residual symbols: Circle for pressure, square for 
x-momentum, triangle for ^/-momentum. 
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Figure 33: Convergence behavior of skin-friction drag coefficient for laminar flow over 
Karman-Trefftz airfoil (a = 0°, Re c = 200, OB = 20 chords). 
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